#include "stdafx.h"
#include "common.h"
#include "marcher.h"

double Godunov(std::vector<Cell*>& src,int i,int j, int w, int h, float dir)
{

	double xum,xup,yum,yup;
	double dx=1.0;
	double dy=1.0;


	xum=fabs(src[AT(i-1,j)]->dist);
	xup=fabs(src[AT(i+1,j)]->dist);

	yum=fabs(src[AT(i,j-1)]->dist);
	yup=fabs(src[AT(i,j+1)]->dist);


	unsigned key=0;


	if (src[AT(i-1,j)]->typ==Accepted) key|=0x1;
	if (src[AT(i+1,j)]->typ==Accepted) key|=0x10;

	if (src[AT(i,j-1)]->typ==Accepted) key|=0x100;
	if (src[AT(i,j+1)]->typ==Accepted) key|=0x1000;

	double u;

	switch (key)
	{

	case 0:
		{
			_ASSERTE(false);
			return 1.0;
			break;
		}
	case 1:
		{
			u=xum+dx;
			break;
		}
	case 0x10:
		{
			u=xup+dx;
			break;
		}
	case 0x11:
		{
			u= (dx + xum + xup)/2.0;
			break;
		}
	case 0x100:
		{
			u= yum+dy;
			break;
		}
	case 0x101:
		{
			double d=dx*dx + dy*dy - xum*xum + yum*(2*xum - yum);
			if (d<0) d=0;
			u = (dx*dy*sqrt(d) + dx*dx*yum + dy*dy*xum)/(dx*dx + dy*dy);
			break;
		}
	case 0x110: //yum,xup
		{
			double discr=dx*dx + dy*dy - xup*xup + yum*(2*xup - yum);
			if (discr<0) discr=0;
			u = (dx*dy*sqrt(discr) + dx*dx*yum + dy*dy*xup)/(dx*dx + dy*dy);
			break;
		}
	case 0x111: //yum,xup,xum
		{
			double d=dx*dx + 4*dy*dy - xum*xum + (2*yum - xup)*(2*xum + xup - 2*yum);
			if (d<0) d=0;
			u = (dx*dy*sqrt(d) +
				dx*dx*yum + 2*dy*dy*(xum + xup))/(dx*dx + 4*dy*dy);
			break;
		}
	case 0x1000:
		{
			u= yup+dy;
			break;
		}
	case 0x1001: //yup+xum
		{
			double d=dx*dx + dy*dy - xum*xum + yup*(2*xum - yup);
			if (d<0) d=0;
			u = (dx*dy*sqrt(d) + dx*dx*yup + dy*dy*xum)/(dx*dx + dy*dy);
			break;
		}
	case 0x1010: //yup+xup
		{
			double d=dx*dx + dy*dy - xup*xup + yup*(2*xup - yup);
			if (d<0) d=0;
			u = (dx*dy*sqrt(d) + dx*dx*yup + dy*dy*xup)/(dx*dx + dy*dy);
			break;
		}
	case 0x1011: //yup+xup+xum
		{
			double d=dx*dx + 4*dy*dy - xum*xum + (2*yup - xup)*(2*xum + xup - 2*yup);
			if (d<0) d=0;
			u = (dx*dy*sqrt(d) + dx*dx*yup +
				2*dy*dy*(xum + xup))/(dx*dx + 4*dy*dy);
			break;
		}
	case 0x1100:
		{
			u = (dy + yum + yup)/2.0;
			break;
		}
	case 0x1101: //yum+yup+xum
		{
			double d=4*dx*dx + dy*dy - 4*xum*xum + (yum + yup)*(4*xum - yum - yup);
			if (d<0) d=0;
			u = (dx*dy*sqrt(d) + 2*dx*dx*(yum + yup) +
				dy*dy*xum)/(4*dx*dx + dy*dy);
			break;
		}
	case 0x1110: //yum+yup+xup
		{
			double d=4*dx*dx + dy*dy - 4*xup*xup + (yum + yup)*(4*xup - yum - yup);
			if (d<0) d=0;
			u = (dx*dy*sqrt(d) + 2*dx*dx*(yum + yup) +
				dy*dy*xup)/(4*dx*dx + dy*dy);
			break;
		}
	case 0x1111:
		{
			double discr=dx*dx + dy*dy - xum*xum - (xup - yum - yup)*(2*xum + xup - yum - yup);
			if (discr<0) discr=0;
			u = (dx*dy*sqrt(discr) + dx*dx*(yum + yup) +
				dy*dy*(xum + xup))/(2*(dx*dx + dy*dy));
			break;
		}

	}

	_ASSERTE(u>=0);
	return fabs(u)*dir;
}


